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ABSTRACT 

High-energy radiation processes in compact cosmic objects are often expected to 
have a strongly non-linear behaviour. Such behaviour is shown, for example, by 
electron-positron pair cascades and the time evolution of relativistic proton 
distributions in dense radiation fields. Three independent techniques have been 
developed to simulate these non-linear problems: the kinetic equation approach; the 
phase-space density (PSD) Monte Carlo method; and the large-particle (LP) Monte 
Carlo method. 

In this paper, we present the latest version of the LP method and compare it with 
the other methods. The efficiency of the method in treating geometrically complex 
problems is illustrated by showing results of simulations of ID, 2D and 3D systems. 
The method is shown to be powerful enough to treat non-spherical geometries, 
including such effects as bulk motion of the background plasma, reflection of 
radiation from cold matter, and anisotropic distributions of radiating particles. It can 
therefore be applied to simulate high-energy processes in such astrophysical systems 
as accretion discs with coronae, relativistic jets, pulsar magnetospheres and gamma- 
ray bursts. 

Key words: plasmas - radiation mechanisms: non-thermal - radiative transfer - 
methods: numerical - galaxies: jets - gamma-rays: bursts. 


1 INTRODUCTION 

The spectral and temporal properties of X-ray and gamma- 
ray emission from many cosmic compact objects suggest the 
involvement of non-thermal radiation processes. Non- 
thermal processes give rise to the radiation from relativistic jets 
in active galactic nuclei (AGNs; see, e.g., review by Sikora 
1994), from pulsars (Daugherty & Harding 1982; Chen & 
Ruderman 1993), and, possibly, from gamma-ray bursts 
(Rees & Meszaros 1992; Meszaros & Rees 1993). They may 
also occur in accretion disc coronae in galactic X-ray 
binaries (Lingenfelter & Hua 1991; Sunyaev et al. 1992) and 
in AGNs (Zdziarski, Lightman & Maciotek-Niedzwiecki 
1993; Madejski et al. 1994). The interpretation of observed 
radiation spectra from these sources is not unique, however, 
and the acceleration mechanisms of relativistic particles are 
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unknown. One way to obtain a better understanding of the 
physical nature of these objects is to perform numerical 
simulations of the radiation processes in different scenarios 
and to compare the results with observations. 

Unfortunately, this task is not easy since the non-thermal 
processes in highly compact objects can be very non-linear. 
For example, consider a type of relativistic particle that radi- 
ates due to interactions with photons (e.g. Comptonization 
by electrons or photomeson production by protons). The 
particles can interact with ‘external’ (ambient) photons, as 
well as with ‘internally’ generated photons produced either 
by the particles themselves or by cascade processes initiated 
by the particles. Furthermore, the number density of radiat- 
ing particles can be strongly increased by pair production 
resulting from absorption of gamma-rays by X-rays, with 
these photons themselves being produced by relativistic 
pairs. 

Non-linear radiation processes have been simulated suc- 
cessfully by solving a system of kinetic equations (Fabian et 
al. 1986; Kirk & Mastichiadis 1989; Coppi 1992; Coppi, 
Blandford & Rees 1993). The kinetic equation approach 
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readily achieves good particle statistics and is very efficient 
for geometrically simple models. However, the efficiency of 
the method falls dramatically with the increasing dimension- 
ality of the system (Xu 1994), because it requires discretiza- 
tion of phase space as well as computations of interactions 
between each pair of momentum cells. In addition, the end 
products of the interactions must be distributed among many 
momentum cells within the allowed kinematic limits. 

Monte Carlo methods can be much more economical as 
regards the number of computational steps, since it is not 
necessary to calculate all possible interactions among par- 
ticles in different regions of phase space. Instead, one selects 
and follows specific particles in order to characterize the 
overall distributions. The trick is then to find the best way to 
represent the particle distributions statistically. To date, two 
approaches have been used to accomplish this: the phase- 
space density (PSD) array representation and the large- 
particle (LP) representation. 

In the PSD approach, the system is represented in the 
same way as in the kinetic equation method - by densities of 
interacting particles in cells of a discretized phase space. The 
transition of a particle between spatial cells is simulated using 
an escape probability formalism, and particle transitions in 
momentum space due to interactions arc simulated using 
standard Monte Carlo methods. So far, this approach has 
been used only to simulate homogeneous spherically sym- 
metric systems (Stern 1985; Xu 1994). 

In the LP approach, the system is represented by an array 
of large particles’, each of which corresponds to a group of 
real particles sharing the same physical parameters - particle 
type, position and momentum vector. Compared to the PSD 
approach, the LP method uses a more flexible representation 
of the particle density in phase space. Associated with each 
LP is a statistical weight, which is proportional to the number 
of real particles represented by the LP. Adjustment of the 
weight allows one to improve the statistical representation in 
those regions of phase space that contain relatively fewer 
particles. Such a representation can be especially useful in 
non-thermally radiating systems, where sparsely populated 
regions of phase space may contain most of the energy. This 
type of representation was used to model non-linear systems 
with Coulomb-type interactions, such as gravitating systems 
or charged particle beams in accelerators. To treat non-linear 
radiation processes near cosmic compact objects using this 
approach, it has been necessary to invent new techniques. 
Such techniques were developed by Stern (1988) for pure 
electromagnetic cascades, and by Stern, Svensson & Sikora 
(1990) and Stern, Sikora & Svensson (1992) for 
hadronic-electromagnetic cascades. 

In this paper we describe the application of the LP method 
to non-linear pair cascades. The LP method has several 
important advantages over other methods: (i) particle trans- 
port is treated naturally, without the use of escape probabili- 
ties; (ii) the efficiency of the code depends very weakly on the 
dimensionality and geometric complexity of the problem; 
and (iii) the representation given by large-particle arrays does 
not depend on model features such as the geometry, the mag- 
netic field structure or the dynamical situation. A code based 
on this representation is easily programmed. The LP method 
is described and compared with other methods in Section 2. 
In Section 3, we perform simulations of non-linear radiation 
processes, both to test the code by comparing with earlier 


work and to show the power of the LP method in solving 
problems never treated before. Finally, in Section 4, we dis- 
cuss and summarize our results. 

2 THE LARGE-PARTICLE METHOD 

Before introducing the LP method it is worth explaining the 
difference between linear and non-linear Monte Carlo simu- 
lations. In a linear Monte Carlo method the objects of a 
simulation are test particles interacting only with a back- 
ground medium, not with each other. The parameters of the 
medium are fixed or updated gradually, taking into account 
the cumulative effects of the incident particles. Since the test 
particle histories arc considered to be independent of one 
another, there is no reason to store information about all par- 
ticle trajectories at the same moment. A well-known applica- 
tion of a linear Monte Carlo approach in astrophysics is the 
work of Pozdnyakov, Sobol & Sunyaev (1977), which deals 
with photon Comptonization. In non-linear situations such as 
those described in the Introduction, however, it is necessary 
to treat the interactions between simulated particles. This 
means that information about all interacting particles must 
be stored simultaneously. We will argue in this paper that the 
LP method is a very efficient technique for organizing and 
storing this information. 

In the LP method, the non-linear system of interacting 
particles is represented by an array of large particles (LPs). 
Each LP is tagged by the type, energy, location in space and 
direction of motion of the physical particles it represents, 
plus a statistical weight cu, which is proportional to the 
number of physical particles represented by the LP. Like the 
physical particles, LPs successively interact with one another, 
changing their attributes and moving in space between inter- 
actions. Because the evolution of all particles must be syn- 
chronized in time, time-steps for the simulation must be 
sufficiently short that changes in the particle populations 
during one time-step are small. We now discuss some import- 
ant techniques used in the LP method: operations with statis- 
tical weights; creation and destruction of LPs; calculation of 
the mean free time of an LP between interactions; choice of a 
target LP; and modifications of the method in cases where 
the straightforward LP technique is inefficient. 

2.1 Statistical weights 

The number of particles in the low-energy range is typically 
several orders of magnitude larger than in the high-energy 
range. Therefore LPs representing particles of different 
energies must have different statistical weights. The LP 
energy weight function for particles of type i is defined to be 

E,(e) = <y f .(f)c, (1) 

where w i is the statistical weight and e is the energy of the 
physical particle represented by this LP. Below we refer to e 
as the LP energy, whereas E t represents the total energy of all 
particles represented by the LP. 

In general, E t can be an arbitrary function of energy, dif- 
ferent for different kinds of particles. The simplest approach 
is to set E i (e) = constant, which means that all LPs share the 
energy of the system equally. The advantages of choosing 
such an energy weight function are that the parts of the par- 
ticle distribution containing most of the energy are best 
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represented statistically, and the number of LPs is conserved 
for a system with fixed total energy. 

An energy-dependent Efe) can, however, be useful; for 
example, if the simulated system includes a luminous and 
stable feature like a blackbody UV ‘bump’ which dominates 
over other parts of the spectrum (as in AGN spectra). In this 
case, we can reduce the number of blackbody’ LPs by using 
a larger E t in the spectral region around the blackbody maxi- 
mum. To cite another example, suppose the spectral radiation 
flux in the system, /* , , is approximately a power-law function 
with a spectral flux index a (i.c. F r <* e~ a ). Here, the choice of 
E i °c e ] ~ a ensures uniform energy resolution over the spec- 
trum. 

2.2 Creation and destruction of LPs 

Since the number of physical particles represented by an LP 
depends on particle energy e (and possibly on time), the 
statistical weight of an LP can change with time, and one has 
to be careful that particle numbers are conserved appropri- 
ately. Particle number conservation is achieved in a time- 
averaged sense, in terms of the expectation values of the 
particle numbers represented by an LP before and after a 
collision. To accomplish this, an LP is either killed or repro- 
duced in a number of identical copies with some probability, 
following each interaction. 

Consider a large particle LP,, with statistical weight w,, 
which changes energy from £, to e 2 (without changing par- 
ticle type) due to an interaction. One can say that LP, is 
replaced by a new large particle, which we denote by LP 2 . 
According to equation (1), the statistical weight of LP 2 is 
given by 


i.c. in general LP 2 will represent a different number of physi- 
cal particles than LP,. Particle number conservation requires 
that the change of the statistical weight be compensated by a 
change in the number of LPs representing the same group of 
physical particles. The mean number of copies of LP 2 pro- 
duced in the interaction must be 


produced with a probability n 3 , otherwise it is produced in 
w 3 copies on average. 

In practice we will consider mainly two-body interactions, 
in which one particle (e.g. LP,) can be considered the inci- 
dent particle and the other (e.g. LP 4 ) the target particle. In 
order to avoid double counting, w f e take into account only 
interactions for which the incident particle has a lower statis- 
tical weight than the target particle, i.e. oj 1 <(o 4 . Since LP 4 
represents a larger number of physical particles than LP,. 
only a fraction a) ] /co A of LP 4 is ‘used up’ in a given interac- 
tion. To conserve target particle number in a statistical sense. 
LP 4 has a probability 1 - (c oJ(o A ) of surviving an interaction. 

Finally, the energy weights of all LPs have a common 
factor w'hich is variable with time, if the total energy of the 
system is variable. This ensures that the number of LPs does 
not exceed the maximum permitted in the array, while main- 
taining enough LPs to obtain close to optimal statistics. For 
example, if the number of required LPs exceeds some techni- 
cal limit, the statistical weight of all LPs can be multiplied by 
a common factor /> 1, and each LP is then killed with prob- 
ability l - f~ \ 


2.3 Maximum cross-section method 

In order to follow the evolution of a given incident LP, we 
must determine its mean free time for interacting with 
another LP and choose a partner from among possible target 
LPs. This could be done by implementing the following 
procedure. 

( 1 ) Calculate the probabilities for the incident LP to inter- 
act with all other LPs contained in the volume V around it, 
taking into account any energy and angular dependence but 
assuming that the spatial distribution of particles in the 
region is uniform. These probabilities are proportional to the 
partial interaction rates, /*-,= (o,d it (e h e„ 6)/ L, where the 
index 7’ refers to the incident LP, the index 7’ refers to the 
target LP, 6 is the angle of interaction, and d lt = o it v i{ is the 
total cross-section multiplied by the relative velocity. 

(2) Calculate the sum = and sample the mean free 
time using the standard Monte Carlo procedure, 

t=-\n{£)/P„ (5) 


_w, e 2 i ) ,,, 

„, = __ = — - (3) 

W 2 E\ E l {E 2 ) 

which generally is not an integer. Therefore we enforce par- 
ticle conservation in a statistical sense, by killing the LP 2 par- 
ticle with probability 1 - n 2 when n 2 < 1, and by copying it 
(on average) n 2 times w'hen n 2 > L For a constant energy 
weight function, the two cases correspond, respectively, to 
energy loss and energy gain by the physical particles repre- 
sented by LP,. 

If a new particle of type /, represented by LP 3 , is produced 
due to an interaction, its expectation statistical weight, w 3 , 
must satisfy 

„■!!!_£> , 4 ) 

co ? e, Eje-i) 

where the energy of the new particle, r 3 , is determined by the 
kinematics of the interaction. For n } < 1, the new particle is 


where £ is a random number between 0 and 1 . 

(3) Choose the type of interaction, the type of target par- 
ticle, and the actual target particle, using the partial interac- 
tion probabilities. 

(4) Simulate the interaction. 

This procedure, while clear and straightforward, is very 
inefficient. It requires 0{N) computational steps (where jV is 
the total number of LPs) per time-step to simulate the evolu- 
tion of each incident LP. Thus the total system simulation 
time scales as 0(N 2 ). 

The efficiency can be improved dramatically by adopting 
a virtual process technique. The trick is to replace the true 
cross-section for each type of interaction /, with a maxi- 
mum total cross-section (including all possible types of inter- 
action), cr^ x , and to assume that v lt = 2c. In this method, the 
maximum partial interaction rates, 
depend only on w ( and do not depend on the properties of 
the incident LP. We define the LP statistical weight cumula- 
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tive arrays S(j) -Y J k^ij a) kl^y where j denotes the LPs in 
order of decreasing energy. The arrays are calculated once at 
the beginning of each time-step, for each type of particle. 
Now, the mean free time can be sampled from equation (5), 
but with P, being replaced by P i max = Y jl P i {i nax , where 

~ ^ max (/max) *^(/min 1 )]» (*) 

where j m jn and y max correspond to the minimum and maxi- 
mum energies of LPs that can participate in the interaction as 
a target [y mm = 1 if there is no energy threshold, and 5(0) = 0 
has to be put in equation (6)]. After the type of target particle 
and type of interaction have been chosen, the index y of a 
possible target LP is searched for from the condition 

S(j) - 5(y min - 1 ) > {[S(/ m J ~ S(j min -D] 

>5(y-l)-S(y min -l). (7) 

After a tentative target, LP y , has been found and all para- 
meters of the interacting particles are known, it is possible to 
calculate the true cross-section, a reah and the relative velo- 
city, v ir of the interaction. This tentative interaction has a 
probability o Teai v [t /d max of being accepted. If it is accepted, 
the interaction is simulated with a standard Monte Carlo 
technique, and the parameters of the LPs are changed 
accordingly; otherwise, the interaction is cancelled and the 
incident LP proceeds on its trajectory without changing 
parameters. 

The maximum cross-section method can easily be under- 
stood in the following way. Let us introduce a virtual channel 
for an interaction that (i) has a cross-section a vjr1 ; (ii) does not 
change the parameters of the interacting particles; and (iii) 
satisfies (a rea , + a vin )v it = er max 2 c. The total interaction rate is 
then proportional to a max , and, after sampling the mean free 
time with d max , we find that the probability of choosing the 
real channel is o ltaX v (l ld max . This is the probability of a real 
interaction mentioned above. 

Now, let us evaluate the simulation efficiency using this 
method. The number of LPs, N, that we deal with should be 
large enough to avoid large statistical fluctuations. The 
number of computational steps for pre-calculating the cumu- 
lative arrays scales as O(N). Since the procedure of finding a 
target LP from equation (7) requires <9[log( N)\ steps for each 
time-step, the total time required for the simulation scales as 
OjNlog(AO). The reduction in the number of computational 
steps from 0( N 2 ) to 0[/V log( N)] improves the efficiency by a 
factor - 10 3 for N ~ 1 .6 x 1 0 4 , the number of LPs used in our 
calculations. 

We note, however, that this method requires extra com- 
putations for virtual-channel interactions. The number of 
attempts to sample an interaction exceeds the number of 
actually simulated interactions by a factor - d max /d real . If the 
behaviour of a rea , is singular (e.g. if there is a resonance), a 
straightforward application of this method would be ineffi- 
cient. One can avoid this inefficiency by treating a resonant 
feature as a separate process that occurs only within a 
narrow energy range. 

2.4 Simulating a spatially inhomogeneous system 

An inhomogeneous spatial distribution of particles is easily 
simulated using the LP method. However, the variable par- 
ticle density must be taken into account when searching for 
target LPs. In the previous section, we noted that a target LP 


is searched for among the LPs in the volume surrounding the 
incident LP. In this volume, the particle distribution is treated 
as uniform. To reproduce an inhomogeneous distribution, we 
simply divide the simulation volume into a number of spatial 
cells. LPs are allowed to interact within the same cell and the 
mean free time is sampled by using the average densities of 
target particles inside a cell. Note, however, that in the LP 
method the trajectories of LPs are followed and the escape or 
transition to another spatial cell is simulated exactly. There- 
fore any cell shape can be used, as long as the spatial cell 
structure is adjusted so that the desired density gradients of 
target particles are obtained. For time-dependent problems, 
spatial cells can be adaptive, i.e. their coordinates and shapes 
can be changed at each time-step. 

In this method, the cumulative arrays are pre-computed 
separately for each cell. This requires little additional 
memory if properly organized. Discretization of the total 
simulation volume into spatial cells hardly affects the 
computing time if the total number of LPs is kept constant. 

2.5 Special modifications of the LP method 

There are situations in which the standard LP method is 
inefficient, due to overly frequent interactions or to poor 
statistics of target LPs. To handle these situations, it was 
necessary to develop a few specialized techniques. Since 
these lie somewhat outside the main line of argument in this 
paper, they are described in the Appendix. 

2.6 Comparison of different methods 

As Xu (1994) argued, the kinetic equation approach is the 
most efficient method for treating 0D (i.e. one spatial zone) 
non-linear problems. To solve the kinetic equations describ- 
ing the time evolution of the system requires ()(N 2 N k Nj 
calculation steps per iteration (i.e. per time-step), where N p is 
the number of bins in momentum space, N k is the average 
number of momentum cells within the kinematic limits for an 
end product of an interaction, and N s is the number of space 
cells. For a typical 0D pair cascade calculation (Coppi 1992; 
Xu 1994), A s = 1, jY p ~ 1 00, N k ~ 10 (the latter estimate is 
very approximate), and the number of calculation steps for 
one iteration of the whole system is - 10\ 

The LP method requires 0[N u Aog{N x P )j calculation steps 
per iteration, where N Li > is total number of LPs. In practice, 
one can neglect the <9(log /V u >) factor, as operations associ- 
ated with this factor are very fast. This estimate does not 
depend on the number of dimensions of the system and satis- 
factory statistics are obtained using only - 10 5 LPs, provided 
that the system is not too complicated geometrically. Thus 
for 0D problems the number of steps is of the same order as 
for the kinetic equation approach. However, each step takes 
longer to compute using the Monte Carlo method, making 
the kinetic equation approach more efficient, in agreement 
with Xu’s conclusion. 

The situation changes drastically where spatial structure is 
taken into account. In a 1 D problem, for example, the kinetic 
equation approach requires at least 1 0 angular bins, making 
/V p £ 10 3 . The required number of spatial cells is at least a 
few. Thus the number of calculation steps for one iteration is 
- 1() 7 , or even more, since N k can grow when angular bins 
are introduced. This is still a reasonable value for simulations 
done on PC-486 computers. However, the number of steps 
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required by the LP Monte Carlo method is now two orders 
of magnitude smaller. Therefore, provided that the LP com- 
putation time for one step is not longer than 100 times the 
kinetic computation time per step, we conclude that the LP 
method is probably faster already for ID geometries. Note, 
however, that it is not easy to compare the computation time 
per step for the two methods, as interactions are simulated in 
completely different ways. 

We now compare the efficiencies of the two Monte Carlo 
methods. The PSD approach can be characterized by 
()[N p N s N T \og{N p )\ steps per iteration, where N r is the 
number of times the calculation is repeated for each phase- 
space cell. The value of N r is optional, but a larger value will 
yield better statistics (Xu 1994). Since a simulation in the 
PSD representation can be interpreted as the evolution of 
/V jV s N r log(/V n ) LPs, and since each LP requires similar 
computation time in both methods, the relative efficiency of 
the two methods can be deduced simply by comparing the 
number N p Nj\ r with the number of LPs used in the LP 
simulation. [For reasons given above we neglect the log (N p ) 
factor.] 

Suppose that the system is represented in both approaches 
in such a way that N u ,= N p N s N r , i.e. the computation time in 
both representations is of the same order. Which approach 
will then give better statistical and systematic accuracy? The 
answer depends on the kind of information being sought. To 
obtain the best statistics in the parts of the particle distri- 
butions containing most of the energy, the best strategy is to 
use an LP method in which all LPs have the same energy, i.e. 
E i (e) = constant. This automatically introduces a greater 
number of LPs into the regions of phase space containing 
most of the energy. If one is interested in the low-luminosity 
parts of the spectrum, the PSD representation has an advant- 
age due to its use of ‘adaptable’ statistical weights (Xu 1994), 
which provides equal statistical accuracy in equal volumes of 
phase space. The efficiency of simulating the high-luminosity 
regions is lower, however, since most time is spent simulating 
the relatively small amount of energy in the low-density 
regions of phase space. The same goal can be achieved in the 
LP approach by using an energy-dependent Eje) (see Sec- 
tion 2.1), but this is less convenient since it requires a priori 
assumptions about the particle energy distributions. 

Another advantage of the PSD approach is that it uses 
information about the intermediate states of particles which 
evolve rapidly in energy space due to numerous interactions 
(see Xu 1994), and therefore gives a better statistical repre- 
sentation of the particle distribution. This is important if the 
rapidly evolving particles play the role of targets in some 
process, e.g. electrons in the synchrotron self-absorption 
process. Simulation of synchrotron self- absorption in the LP 
approach requires a special treatment, incorporating some 
features of the PSD representation (see the Appendix). 

Summarizing the above, we find that both Monte Carlo 
methods have comparable efficiencies from a statistical point 
of view for problems of low dimensionality. The LP repre- 
sentation is more efficient if we are interested in the high- 
luminosity parts of the spectrum, while the PSD approach is 
more convenient for reproducing the low-luminosity spectral 
regions and to simulate processes with very large cross- 
sections, such as synchrotron self-absorption. 

Note, however, that the equation of N u , and N p N s N r is 
possible only for problems with dimension less than two. For 
2D pair-cascade problems the number of required PSD 


phase-space cells can be estimated to be > 10\ i.e. 100 
energy bins times (20 x 10) angular bins times more than 10 
spatial cells. Note that two angles are necessary even in 2D 
problems. If each cell in phase space is simulated N r > a few 
times, the PSD approach is equivalent in computational 
speed to an LP simulation with a few million LPs. In the LP 
scheme, on the other hand, £ 10 s LPs is sufficient, except for 
problems with very complicated geometry. In addition, PSD 
calculations of higher dimensionality are less efficient 
because a larger fraction of time is spent on simulating 
regions of phase space that are sparsely populated. In order 
to overcome this problem, it is necessary to use complicated 
adaptable mesh techniques, whereas in the LP representation 
an optimal distribution of interactions throughout phase 
space is achieved automatically. 

We now comment on Xu’s (1994) claim that the LP 
method is very inefficient for simulations of dynamical (i.e. 
time-dependent) systems. He argues that the time resolution 
of the LP method is constrained by memory limitations on 
the number of LPs, while statistical precision can be achieved 
in the PSD aproach by performing a large number, /V r , of 
simulations of each phase-space cell. However, consideration 
of the memory requirements does not bear this out. Each LP 
in a 3D problem requires approximately 40 bytes of memory. 
Thus, for 10 5 LPs (in the present paper, 16 384 LPs were 
used), 3,5 Mb of memory is needed and computers like a PC- 
486 can be used. With a few tens of Mb of memory, which 
Xu claimed is necessary for a 2D problem using the PSD 
method, one could simulate 1() 6 LPs. Noting that the escape 
rate of LPs is ~ N lv per light crossing time, R/c , and that 
reasonable statistics are achievable for ~ 10 s LPs, we find 
that a memory capacity of 10 fi LPs seems sufficient to simul- 
ate a system with very good statistical accuracy over a time as 
short as the light crossing time. Thus we claim that, at pres- 
ent, computational speed rather than memory requirement is 
the bottle-neck when simulating non-linear radiating sys- 
tems. 

Let us next consider systematic errors. Both approaches 
have errors associated with the finite value of the time-step, 
A t. This does not present a serious problem, however, 
because these errors ~ O(Az) and can be reduced to any 
small value by choosing a sufficiently small A/. The PSD 
approach involves some additional errors due to the discreti- 
zation of momentum space as well as real space. The effects 
of the former errors can be moderate, but effects of the latter 
can be very unpleasant. A particle in the PSD approach has 
no continuous trajectory. Motion between space cells is 
described by Poissonian statistics, i.e. by a constant prob- 
ability per unit time, dpjdt~ p/r, where the escape time r is 
defined by the size of the space cell divided by the particle 
velocity. As a result, the flight times of a particle through a 
space cell will be distributed with a probability p ~ e~' /r , and 
after crossing many cells the flight times will have a Gaussian 
distribution. This effect is especially troublesome in dynami- 
cal problems. For example, a front of photons emitted by a 
flare will be artificially dispersed in space. The effect will 
cause systematic errors even for steady-state problems - e.g. 
the dispersion of flight times of a photon through some 
region will be overestimated. 

Finally, we emphasize that the space cells that are used in 
the LP approach do not discretize space. Rather, the LP 
method introduces discrete levels of the target density, since 
the targets are assumed to be uniformly distributed within a 
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given cell. Therefore the requirements regarding the size and 
shape of space cells can be considerably weaker than for the 
PSD approach. The main surviving systematic errors are 
those associated with gradients of energy and angular distri- 
butions across a single space cell. Since the locations of all 
LPs are exactly known, however, it is possible to reduce these 
errors by calculating the gradients of target densities within 
each space cell and taking these into account when sampling 
free path time. 

3 SIMULATIONS OF PAIR CASCADES USING 
THE LP METHOD 

3.1 Pair cascades 

The simplest model of a non-thermal pair cascade (see, e.g., 
review by Svensson 1994 and references therein) comprises 
a uniform region of size R , into which monoenergetic pairs 
or electrons with Lorentz factor y max are injected with a 
luminosity L e , along with a luminosity L uv of blackbody 
radiation at temperature T uv (typically in the UV range). 
Since the electrons lose virtually all of their energy by 
Compton scattering the soft photons, the entire luminosity 
L e is converted into X-rays and gamma-rays. 
Electron-positron pairs are typically produced in 
photon-photon interactions between gamma-rays with ener- 
gies e> 1 (in m e c 2 units) and X-rays slightly above the pair 
production threshold l/e. The radiation from the pairs may 
produce further generations of pairs, resulting in a 
pair-photon cascade. 

The self-consistent pair and photon distributions, taking 
into account reprocessing by the pair cascade, depend 
on only four parameters: (i) the injection compactness = 

( LJR)(o T /m c c } ); (ii) the soft (or UV) compactness / l;v = 
(L uv //?)(a x /m c c 3 ); (iii) the soft blackbody temperature kT v fV ; 
and (iv) the injected Lorentz factor y max of the electrons. 
There are essentially three regimes for these pair-photon 
cascades: (i) no pair cascading for t c <a few; (ii) unsaturatcd 
pair cascades for a few ^ / c £ 30; and (iii) saturated pair cas- 
cades for / c >30. For saturated pair cascades, essentially 
every gamma-ray is absorbed in producing an electron-posi- 
tron pair. 

The pairs cool to non-relativistic energies where they 
thermalize and eventually annihilate at the equilibrium 
temperature T c . This temperature is set by the balance 
between Compton heating and inverse Compton cooling on 
the self-consistently determined radiation field. The 
Thomson scattering depth of the cool pairs, r r , becomes 
larger than unity for injection compactness £ 20. Then the 
interactions between the pairs and the photons will have a 
noticeable effect on the radiation field, i.e. the UV and X-ray 
radiation will be Comptonized. The UV blackbody photons 
will be upscattered to produce a steep soft X-ray power-law 
spectrum, while the hard X-rays at energies > l/r 2 r will lose 
energy by scattering on the cool pairs. The latter process will 
cause a spectral break in which the spectrum steepens 
towards higher energies. An important diagnostic parameter 
is the pair yield, PY, which is the fraction of injected power 
that is converted into pair rest mass, Guilbert, Fabian & Rees 
(1983) showed that r T depends on the pair yield according 
to r T ~(/ c PY ) I/2 . The pair yield is an increasing function of 
<? e , but reaches a constant value of about 10 per cent for 
saturated pair cascades (Svensson 1986). 


If the electrons are injected with a power-law distribution 
in energy, the model acquires two additional parameters, the 
slope of the injected distribution, T, and the minimum 
injected Lorentz factor, y mm . Another two parameters enter 
if the magnetic field is taken to be non-zero; these can be 
chosen to be a ‘magnetic compactness’, l B = (B 2 /&n)o T R/ 
m e c 2 , and the source size R (or the magnetic field B). 

3.2 Main features of the LP cascade code 

We have implemented the LP method described in Section 2 
in a code aimed mainly at simulating non-linear electromag- 
netic and hadron-electromagnetic cascades. This code can 
also be used efficiently to study simpler linear phenomena, 
e.g. thermal Comptonization. The prototype of the code, pre- 
sented in Stern (1988), included electromagnetic processes 
only; later versions have been extended to include processes 
resulting from injection of relativistic protons (Stern, Svens- 
son & Sikora 1990; Stern, Sikora & Svensson 1992). The 
types of particles represented by LPs in the code are 
photons, electrons, positrons, protons, neutrons, and 4 Hc 
nuclei and their fragments. The interactions taken into 
account in the present version of the code are Coulomb 
scattering of electrons and positrons, synchrotron radiation, 
synchrotron self-absorption, Compton scattering, 
photon-photon pair production, pair annihilation, pair pro- 
duction in proton-photon interactions, photomeson produc- 
tion (on both protons and neutrons, including p-n charge 
exchange), inelastic proton-proton interactions, and photo- 
disintegration of 4 He and its products (Sikora & Begelman 
1992). In the present version of our code, the pions and 
muons are assumed to be too short-lived to suffer collisions 
before decaying and, therefore, we treat their decay products 
as being produced in situ. This assumption, however, is cor- 
rect only for AGNs and, for simulations of radiation pro- 
cesses in galactic X-ray sources, the collisions of charged 
pions and muons with other particles have to be followed as 
well. 

The microphysics of most electromagnetic processes is 
represented by the exact expressions of quantum electro- 
dynamics (e.g. Jauch & Rohrlich 1976). The only approxi- 
mations we use are the following: the ultrarelativistic 
expressions for synchrotron emission and self-absorption are 
also applied in the semirelativistic regime, and an approxi- 
mate expression for the electron-positron Coulomb scattering 
cross-section is used. The latter assumes that the target par- 
ticle is non-relativistic and that the energy exchanged is negli- 
gible in comparison with the energy of the incident particle. 
It can be shown that these approximations introduce negli- 
gible effects to the results of simulations for typical astro- 
physical situations. Non-electromagnetic processes, such as 
photomeson production, photodisintegration and nuclear 
collisions, presently do not have good theoretical descrip- 
tions and are therefore simulated using experimental cross- 
sections and inelasticities. We do not give references for that 
data here because the main purpose of this paper is the 
demonstration of the efficiency of the LP method for simple 
scenarios involving only electromagnetic processes. 

Simulations of all processes, including the spiralling of 
charged particles in large-scale magnetic fields, are imple- 
mented using full 3D kinematics. Charged particles, being 
coupled to the background plasma through chaotic magnetic 



Non-linear high-energy processes near compact objects 297 


fields, also participate in externally imposed bulk motions. 
Multiple spatial cells are used to simulate inhomogeneities in 
the target medium. A few easily reprogrammable subroutines 
constitute the model-dependent part of the code, while the 
main part of the code is universal. 

The efficiency of the LP code can be characterized as fol- 
lows. For 2 14 ( = 1 6 384) LPs the CPU time for simulating the 
evolution of a pair cascade over the time interval ~ Rjc, 
where R is the characteristic size of the system, varies 
between 0.1 and 0.5 h (depending on the compactness para- 
meter) on computers such as a Sun-4 or a 486DX 2-66 PC. 
Usually, it is necessary to follow the system evolution for 
~ 10-20 (R/c), in order first to reach a steady state and then 
to achieve a statistical quality similar to that shown below in 
the figures. 


3.3 Comparison of LP and kinetic equation results 

Here wc compare the results of our LP cascade code with the 
kinetic equation results of Coppi (1992). There are earlier 
works on pair cascades carried out within the framework of 
the kinetic equation approach (e.g. Fabian et al. 1986; Light- 
man & Zdziarski 1987; Svensson 1987; Done & Fabian 
1989; Ghisellini 1989), but Coppi (1992) contains the most 
detailed treatment of the microphysics and covers a larger 
variety of cascade parameters. 

Physical space in our calculations is represented by a 
homogeneous sphere. This corresponds most closely to the 


one-zone approximation, i.e. a uniform medium with an 
escape probability, that was used by Coppi (1992). In this 
series of calculations we follow Coppi in assuming external 
injection of pairs. The contribution of the injected pairs to 
the total pair yield is small in all cases except that repre- 
sented in Fig. 2, where the intense annihilation line is an arte- 
fact of the injection assumed by both Coppi and ourselves. 

3.3. 1 Simple Compton pair cascades 

High-energy processes near accreting compact objects often 
occur in a background of thermal radiation from optically 
thick matter. Since this thermal radiation emerges in the UV 
band in AGNs, most previous work on pair-cascade models, 
including that of Coppi (1992), has included such a UV radi- 
ation component. In these models, the main source of pair- 
producing photons is Comptonization of the UV photons by 
relativistic electrons. 

We compare our results for Compton cascades with those 
of Coppi, for three cases. Three of the input parameters are 
kept constant - monoenergetic pair injection w ith y max = 10\ 
/c7jj V = 10“ s m c c' 2 and <? llv // e = 2.5 - while takes the 
values 10, 100 and 1000. The escaping radiation spectra for 
the three values of <? e arc presented in Fig. 1 (solid curves), 
together with the spectra from fig. 1 in Coppi (1992; dotted 
curves). Three derived parameters, the total scattering opti- 
cal depth r T , the pair yield and the pair temperature T e , are 
compared in Table 1 with the corresponding results from 
table l in Coppi ( 1992). 



Log(e r /m„c 2 ) 

Figure 1 . Spectra of escaping radiation from LP calculations (solid lines) and from kinetic equation models (Coppi 1992; dotted lines), for 
unmagnetized pair cascades with compactnesses / c =10, 100, 1000. /,, v // c = 2.5 for all cases. The sharper-than-exponential slope of the 
blackbody spectrum visible in our results is due to our artificially cutting off the exponential tail of the Planck distribution. The blackbody 
photon energy distribution has been corrected to the full Planck spectrum in calculations presented in Figs 8- 1 2. 
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Table 1 . Comparison of LP and kinetic equation simulations. 


l e = 1000 

r T 

Pair Yield 

kT t /m e c 2 

LP (this paper) 

13.2 

0.121 

4.5 x 10“ 

kin. (Coppi 1992) 12.8 

0.12 

6.2 x 10" 

U = ioo 

LP 

3.41 

0.082 

2.3 x 10" : 

kin . 

3.34 

0.087 

2.2 x 10’ : 

o 

II 

u 

LP 

0.405 

0.014 

4.9 x 10~ ; 

kin. 

0.502 

0.023 

5.7 x 10" ; 


The differences between the two sets of results, while 
relatively small, are physically significant, in the regime of 
unsaturated pair cascades near / e ~ 10, the results are known 
to be sensitive to details of the model, such as how the 
opacities are calculated and how the radiative transfer is 
treated. The differences in the radiation spectra for / c = 10 
may not seem large. A small difference in the X-ray spectrum 
(and thus the optical depth to pair production), however, 
causes a small difference in the gamma-ray spectrum, which 
in turn leads to a substantial difference in the number of 
pairs produced and thus in the pair yield. Thus Table 1 
shows a difference in the pair yield of - 60 per cent. 

The main cause for these differences lies in the different 
treatments of photon escape used in the two methods. Coppi 
(1992) uses a simple escape probability for his one-zone 
model. He implements the prescription for the mean escape 
probability per second, P csc (i.e. the inverse of the mean 
escape time), originally introduced by Lightman & Zdziarski 
(1987): 

Asc* -«(c)/fl(*) = [l + t T /(f)] _1 , (8) 

where f(e) = 1 for e ^0.1, f(e)-( 1 - £)/0.9 for 0.1 <e< 1, 
and /(e) = 0 for £ > 1. In the limits of r T ^> 1 and <3C 1, P csc is 
known analytically for various geometries and spatial distri- 
butions of photon sources and scatterers. The choice made 
by Coppi corresponds to a homogeneous source with slab 
geometry. For the spherical homogeneous source we con- 
sider, P csc = 4c/3 R for r T <5C 1 , and 5c/t t R for r T » 1 and 
£<C1 (Sunyaev & Titarchuk 1980). Fitting the LP simula- 
tions, we obtain 

Asc oc (3/4 + 0.188t x )- 1 (9) 

for e <5C 1, in good agreement with the analytic results. 

Coppi’s use of equation (8) underestimates P csc for a 
spherical source, resulting in a larger X-ray density which 
implies, in turn, a larger gamma-ray opacity and the turn-on 
of pair production at a lower / e . Thus, the turnover of the 
photon spectrum for / e =10 in Fig. 1 occurs at a lower 
energy, and both r T and the pair yield (in Table 1) are larger 
in Coppi’s simulation. 

For larger values of the compactness, the cascade is satu- 
rated. Then r T and the pair yield are insensitive to the details 
and all simulations should give the same results, independent 
of method. This is clearly seen in Table 1 for / e = 100 and 
1000. For a fully saturated case, there is effective downscat- 


tering and essentially all of the reprocessed gamma-ray 
power emerges at £ £ 1 /r x . The radiation spectra obtained 
using different simulation methods should agree at these 
energies, simply due to energy conservation. At larger 
photon energies, photons escape from a surface layer of unit 
optical depth and there may be large differences in results 
obtained using different prescriptions for P csc . Because 
Coppi underestimates P esc for £ < 1 , photons spend a longer 
time in the source, allowing a larger fraction of the photons 
in this energy range to downscatter to £ ~ l/r£ before escap- 
ing. Thus Coppi’s spectrum steepens by unity in the spectral 
index at the hard X-ray break at £ ~ 10 2 for 10\ while 
our corresponding spectrum steepens by 1 /2, as is expected 
for the case of a uniform distribution of photon sources and 
scatterers (Sunyaev & Titarchuk 1980). At £> 1, the opacity 
due to photon-photon pair production dominates at large / c . 
To have too small an X-ray density for l/r|<£< 1 causes 
the opacity for gamma-rays to be smaller, and the flux of 
photons escaping with 1 <e<t\ will be too large. This is 
seen in Fig. 1 for / c = 1 0\ where Coppi's flux slightly exceeds 
ours at £> 1. 

The temperature of the cooled pairs given in Table 1 is 
determined from a balance mainly between Compton cool- 
ing and heating, both of which are sensitive to the spectrum. 
As the spectra obtained with the two methods differ, so will 
the temperatures. Finally, we remark on an artificial feature 
visible in Figs 1-6. The sharper-than-exponential slope of 
the blackbody spectrum is not physical but results from our 
using a sharp cut-off at the high-energy end of the Planck 
distribution. This crude approximation is corrected in the 
models calculated in Sections 3.4.5 and 3.4.6, and in Figs 
8-12. 

3.3.2 Models including synchrotron radiation 

We have modelled the case presented by Coppi (1992) in 
which the electrons also lose energy through synchrotron 
emission, with the following parameters: magnetic field 
/? = 3()0G, magnetic compactness £ H = 2.9 {R= 10 ls cm), 
UV compactness / uv = 6, UV temperature kT vv = 1.07 x 
l()- 5 wi c c 2 , power-law injection spectrum with slope T = 2.4, 
maximum electron energy = 1 0\ and different minimum 
injection energies. The agreement is good (Fig. 2), except 
that the ratio of Comptonization to synchrotron emission is a 
bit lower in our case, the reason for this being the higher 
escape probability in our models. 
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Figure 2. Spectra of escaping radiation from magnetized pair cascades, calculated by the LP method (solid lines) and the kinetic equation 
method (Coppi 1992; dotted lines). For all models. fl = 300G, /„ = 2.9 (R= 10 15 cm), i c = 30, / uv = 6, and we assume power-law injection of 
relativistic electrons with slope T = 2.4, y m „= 1 000, for different values of y min . Curves are normalized by factors shown in parentheses. 


3.4 New results from LP calculations 

In this section we present a series of examples illustrating the 
capabilities of the LP method. We successively model 
inhomogeneous pair cascades, dynamical effects and com- 
plex geometries. 

3.4. 1 Multishell representation for uniform injection 

We start with the same case of a pure Compton cascade and 
uniform injection as in Section 3.3.1, but with the source 
volume divided into seven concentric shells, five internal 
shells where relativistic electrons are injected uniformly and 
two external shells with zero injection. The internal shell 
stratification gives a better representation of the spatial dis- 
tribution of pairs and the spatial and angular distributions of 
target photons. This allows us to test whether these distri- 
butions are roughly uniform inside the source, which might 
be expected in the case of uniform injection. To emphasize 
the effect of pair creation in the space surrounding the injec- 
tion region, we assume a very high electron injection energy, 
y m;ix = 1() 6 , with / e = 100 and / uv =400. In Fig. 3 we com- 
pare the results with the equivalent one-shell model, in which 
energy is injected within the same volume as for the multi- 
shell model. The difference is visible only in the highest- 
energy part of the spectrum, where the multishell model 
gives an exponential cut-off, whereas the one-shell model 
predicts a power-law tail. The reason for the difference is 
that the multishell model takes into account pair creation 
outside of the injection region; due to the very high electron 
injection energy, gamma rays are effectively absorbed in this 
region. The X-ray part of the spectrum is insensitive to 
spatial stratification, which can be explained by the fact that 
the pair distribution inside the injection region does turn out 
to be almost uniform. We can therefore conclude that a one- 
shell approximation for a region of uniform injection is 
adequate. 


3. 4 . 2 Non - un iform injection 

We next consider a case in which relativistic electrons 
and UV radiation are injected non-uniformly. We assume 
that both injection functions are distributed according to 
dL/d/?~ 1 /R 2 and that the ratio of the maximum injection 
radius to the minimum injection radius is 5. The escaping 
spectrum is presented in Fig. 4, together with the spectrum 
for uniform injection (as in Section 3.4.1), both calculated 
assuming the same multishell structure. The spectra are 
almost identical in the photon energy range e< 1, implying 
that the uniform injection approximation is probably 
adequate for making predictions about X-ray spectra. 

At higher energies, the non-uniform injection model gives 
a flux that is lower by a factor - 2. This can be explained in 
terms of gamma-ray absorption. In the non-uniform injection 
case, gamma-rays are weighted towards the central shells. 
Since the compactness is higher there, the optical thickness 
for gamma-ray absorption by pair production is also higher 
and the net flux of escaping gamma-rays is smaller. 

3.4.3 Additional heating of thermal pairs 

Motivated by results from the OSSE instrument on Compton 
Gamma-Ray Observatory (Johnson et al. 1994), showing that 
spectra of Seyfert galaxies often have a sharp break at ener- 
gies corresponding to temperatures much higher than the 
expected Compton temperatures in pair-cascade models, we 
calculated a model in which a large fraction of the energy is 
injected by direct heating of the thermal pairs (i.e. those that 
have cooled from relativistic energies but have not yet annihi- 
lated). Such heating was considered by Zdziarski et al. 
(1993). It can be caused by Coulomb interactions with semi- 
relativistic protons or just by an acceleration process of finite 
efficiency, as was suggested by Zdziarski et al. (1993). High 
temperatures of scattering electrons can also be mimicked by 
high velocities of turbulent elements, provided that the 
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Figure 3. The effects of using multiple spatial shells on pair-cascade spectra calculated using the LP method, for a high electron injection 
energy, y max = 10\ with <f t , = 100, / uv = 400. The solid curve is for the one-shell model and the dotted curve is the equivalent model using seven 
shells. 



Figure 4. 1 he effects of non-uniform relativistic electron injection on escaping radiation spectra are shown for models with seven spatial shells. 
Solid curve - electrons are injected according to the radial distribution dT mj /dr~ r 1 within the limits 0.2 < r< 1. The parameter equals the 
compactness referred to r = I, and L l v =4L v . Dotted curve - model with the same parameters, but with uniform electron injection. 


characteristic size of turbulent cells does not exceed the 
photon mean free path. 

In our calculations we artificially set a constant tempera- 
ture and performed simulations for the same shell structure 
as in the previous two examples. The power that must be 
supplied to the pairs to support the assumed temperature is 
found to be - 3 times higher than the power injected through 
relativistic electrons. The results are presented in Fig. 5 and 
compared there with the purely thermal model. 


The lack of a UV bump means that practically all UV 
photons undergo Comptonization. In order to reproduce 
observed AGN spectra, which generally possess prominent 
UV bumps, one would have to appeal to additional geometric 
complexity in the source. For example, one could posit that 
the UV source is spatially separated from the cascade region, 
or that the active corona only partially covers the accretion 
disc. Partial covering could result from the local character of 
coronal activity (magnetic flares), and the domination of 
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thermal over non-thermal Comptonization could occur if 
most of the energy dissipated in flares goes to heat the 
plasma. 

3.4.4 Effect of radial bulk motion 

To study the possible effects of a radial bulk motion on the 
spectrum, we made calculations for the same set of para- 
meters as in Section 3.4.3, adding a constant inward radial 
velocity (accretion) and assuming an r 2 distribution of 
electron injection. In this case the temperature is fixed in the 


rest frame of the pair flow. Escaping radiation spectra for 
v/c= -0.1 are shown in Fig. 6 in comparison with the same 
case for r = 0. 

As expected, accretion hardens the spectrum. This effect 
is essentially a form of adiabatic heating, in the sense that the 
compressional force does work against the photon gas. In 
terms of radiation transport it can be described as first-order 
Fermi acceleration of photons in the converging flow: a 
photon is more likely to be scattered by an oncoming 
electron than by a receding one, leading to a mean increase 
of photon energy per scattering. The mean number of scat- 
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Figure 5. The effects of additional pair heating. The temperature of the thermal pair distribution was artificially supported at the constant 
value kl '= 0.03 m K c 2 . Solid curve - Comptonized cascade spectrum, for injected power with compactness parameters A • = 100, A v = Id. The 
total escaping luminosity L is three times larger than the injected power, with the additional energy supplied through pair heating. Dotted 
curve - Comptonization of the UV Planck spectrum by heated thermal pairs alone, for A = d, Arv = Id and A7’=0.03w c c 2 . The density of 
thermal pairs was set equal to that obtained in the non-thermal cascade calculation. 



Figure 6. The effect of radial bulk motion with constant inflow velocity v = -0.1 c on the escaping spectrum at fixed pair temperature 
A7' = 0.()3m L .c 2 , for A = 100, A-v = Id and an r 2 radial distribution of electron injection. Dotted curves - no photon absorption in the centre is 
simulated. Solid curves - photons are absorbed at r<0.1 5. 
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terings becomes 26.4, versus 1 1.3 for v = 0 , but the optical 
depth is almost the same as the model with v = 0. The escap- 
ing luminosity is 17.6 times larger than the non-thermal 
power, and a thermal break in the escaping spectrum is 
almost transformed in shape into a Wien peak. 

In the case shown by dotted lines in Fig. 6, all photons are 
allowed to escape. In reality, a large fraction of the photons 
produced in the central region of an accretion flow could be 
absorbed by the black hole. This case is approximated by our 
second model (solid lines in Fig. 6), in which all photons and 
pairs that venture within 0.15/? of the centre are killed, 
where R is the radius of the region where energy injection is 
taking place. For the case of accretion with v= - 0.1, the 
effect is very strong - the escaping luminosity drops almost 
by an order of magnitude and the mean number of scatter- 
ings before escape drops by more than a factor of 2. For 
v = 0, the effect of absorption at the centre is smaller but still 
significant. 


3.4.5 Two-dimensional simulation of a non-thermal corona 
above an accretion disc 

Reflection and reprocessing of X-ray radiation by optically 
thick gas (e.g. an accretion disc) is thought to be an important 
feedback mechanism in the production of high-energy 
spectra in AGNs (Haardt & Maraschi 1993). All geometrical 
configurations that include a disc, except perhaps an infinite 
slab geometry, are at least two-dimensional. We calculated 
two simple models to illustrate how the LP method could 
handle both the higher dimensionality and the feedback self- 
consistently. As in previous models we assume that the ther- 
mal particles in the corona are subject to an additional heat 
source. Instead of fixing the temperature as in the previous 
two subsections, however, we fix the power supplied to the 
thermal electrons L, and calculate the temperature self- 
consistently, as well as the Compton reflection component 
from the disc and the final escaping spectrum. 


As a simplified model for reprocessing, we follow each 
photon entering the disc, taking into account only Compton 
scattering (photoelectric absorption is neglected in these 
calculations, but is included in the calculations presented in 
Section 3.4.6). If the photon is reflected after fewer than 10 
Compton scatterings, we follow its further trajectory, starting 
with the energy at which it leaves the disc. The energy depo- 
sited in the disc is emitted in the form of a blackbody UV 
photon from the point at which the incident photon entered 
the disc. If a photon remains in the disc through 10 or more 
Compton scatterings, we assume that the whole photon 
energy is re-emitted in the form of UV blackbody radiation 
with a specified temperature. 

We consider two geometrical configurations: (a) hemi- 
sphere + disc; and (b) finite cylinder + disc. 

(a) Hemisphere ’ + disc. This model consists of an infinite 
plane parallel slab (‘disc’) surmounted by a finite hemispheri- 
cal region of hot plasma above the disc surface (Fig. 7a). The 
reflection and re-emission from the disc beyond the hemi- 
sphere is taken into account. This configuration can be asso- 
ciated with two physical scenarios - a global hot plasma 
corona (in which case the black hole must be located in the 
centre of the hemisphere) and a local bulge of hot plasma on 
the surface of the disc. The spatial cell configuration is shown 
in Fig. 7. 

We performed two calculations for this geometry. The first 
is a test ol pure thermal Comptonization, with reprocessing 
and reflection. The hemisphere was homogeneously filled by 
electrons with an optical depth r T = 2 across the radial dis- 
tance. We set the ratio T t /TuY = 25, so that the intrinsic 
blackbody emission of the disc is negligible and the main 
share of UV radiation is produced as the result of the feed- 
back. The absolute value of the power does not matter as we 
have no photon-photon interactions and the result can be 
scaled to any total power as well as to any size of the source. 
The resulting spectrum in direction range cos 0 >0.5 
(measured with respect to the disc normal) is shown in Fig. 8 




Figure 7. Geometric arrangements of spatial cells used to model hot plasma + disc systems. The disc is modelled by an infinite plane slab; 
numbers mark the spatial cells in the regions of hot plasma. All sizes are given in terms of the radius of the hot plasma region, I. (a) Hot 
plasma occupies a hemisphere sitting on the disc surface. The coordinates of the LPs are two-dimensional, (b) Hot plasma occupies a cylinder 
of length 4, the lower edge of which is situated a distance h above the disc surface. In this case the coordinates are three-dimensional, but 
simulations are performed for only one-half of the cylinder. Full dimensionality is achieved by mirror reflection of all photons crossing the 
median plane. 
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Figure 8. Escaping photon spectra for the case of hemisphere + disc in the angular range cos #>0.5, where 6 is the angle from the disc 
normal (see text and Fig. 7). Solid curve - combined non-thermal and thermal power supplies: / e = 10, / t = 5(). Dotted curve - only thermal 
power supply to ambient electrons: / t = 50, r , = 2. 


(dotted line). Its slope corresponds to an energy spectral 
index a ~ 1.09. Note that in this geometry it is impossible to 
get much harder (a<l) spectra due to strong feedback, 
unless r T < 1 (Haardt & Maraschi 1993). The electron tem- 
perature in the hemisphere is found to vary as a function of 
position, increasing from kT= 0.06 m c c 2 at the bottom to 
0.14ra c ,c 2 at the top. 

In the second model we do not assume the existence of 
any ambient thermal plasma. Instead, only relativistic 
electrons are injected. The injection parameters in this model 
are / c - 10, /, = 50, / uv = 2 and y c = 10 6 . The resulting spec- 
trum is shown in Fig. 8 (solid line). Here the spectral index is 
a =1.06, and the temperature of the thermal pairs in the 
hemisphere increases from 0.10 at the bottom to 0.17 at the 
top. The pair density distribution is close to uniform, with a 
maximum density contrast (ratio of the maximum density to 
the minimum density) of 1 .34. 

The anisotropy of the radiation can be estimated from Fig. 
9, which shows the spectrum of the latter model for two 
ranges of polar angle, cos 0 > 0.5 (solid line) and cos 0 < 0.5 
(dotted line). 

(b) Cylinder ( magnetic flux tube) + disc. The feedback 
can be weaker if the region of hot plasma is separated by 
some distance from the surface of the disc. We represent this 
kind of geometry, which could be associated with a magnetic 
flux tube in the corona, by a cylinder of radius r whose axis 
lies parallel to the disc surface. The length of the cylinder is 
4 r, and the gap h between the cylinder and the disc is 0.5 r. 
The geometry of the spatial cells is shown in Fig. 7(b). The 
cylinder was not divided into space cells along its axis, imply- 
ing a 2D target density distribution within the cylinder, but 
particle escape was treated fully three-dimensionally. 

As above, we assumed LJL m = 25, where L uv is the 
intrinsic disc radiation penetrating into the cylinder. We cal- 
culate the same variants as in the previous geometry. For the 
case of pure thermal Comptonization when electrons fill the 
cylinder homogeneously with r=2 (from the axis to the sur- 


face), the resulting spectral index is a = 0.83 (Fig. 10, dotted 
line) and the temperature increases from the side closest to 
the disc to the side furthest from the disc, from 0.06 to 0.10. 
The UV bump is dominated by the reprocessing of hot 
plasma radiation by the disc. 

In the second case we use injection parameters = 1 0 and 

= 50 (obtained by dividing the appropriate luminosity by 
the radius of the cylinder, so that the total power and effec- 
tive compactness are higher than in the case of the hemi- 
sphere). The resulting spectral index is a = 0.86. Again, 
temperature increases from the bottom to the top, from 0.08 
to 0.15. The pair distribution is almost uniform, with a maxi- 
mum density contrast of 1.24. 

The calculated feedback coefficient for this geometry is 
L fo/L U)t = 0. 1 4, where L ^ is the luminosity re-emitted by the 
disc and entering the hot plasma region, and L U)t is the total 
power supplied to the hot plasma. This feedback coefficient 
is probably the main parameter defining the maximum X-ray 
slope that can be achieved in a given geometry w ith multiple 
thermal Comptonization. For example, for the case of a spec- 
trum with a <0.8 and with a thermal break such as that of 
NGC 4151, we can conclude that the separation of emitting 
region from the disc is higher than in the case considered 
here. To demonstrate the effect of spatial separation between 
the region of hot plasma and the disc, we simulated a hybrid 
model with the same parameters as in the previous example, 
except for a larger gap h = 1 .5 r between the cylinder and the 
disc. The result is presented in Fig. 11 (dotted line), in com- 
parison with that for h = 0.5 r (solid line). The spectral index 
for the case h = 1.5 r is a = 0.7, and the temperature and pair 
density are ~ 10 per cent higher than in the case h = 0.5 r. 

3.4.6 Thermal pair production and effect of photoelectric 
absorption in the disc 

Thermal pair production is one of the effects which can 
regulate the optical depth of the hot plasma (Ghisellini & 
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Figure 9. Escaping photon spectra for the case of hemisphere + disc, combined non-thermal and thermal power supplies as in Fig. 8, for two 
angular ranges. Solid curve - cos 0> 0.5. Dotted curve - cos 8 < 0.5. 



Log(e r /m e c 2 ) 

Figure 10. Escaping photon spectra tor the case of cylinder 4- disc with =0.5, in the angular range cos 0>O.5 (see text and Fig. 7). Solid 
curve - combined non-thermal and thermal power supplies: / c = 10. /, = 40. Dotted curve - only thermal power supply to ambient electrons: 
A = 50, r, = 2 over the distance of the cylinder's radius. 


Haardt 1994, and references therein). We demonstrate the 
effects of thermal pair production for the same geometry 
(cylinder + disc, h = Q.5r) as in the previous section. The 
thermal power supply corresponds to / t = 50. This is the 
same as the total power in our hybrid model, the difference 
being that the non-thermal component is now absent. No 
pre-existing matter has been assumed to occupy the cylinder, 
so that the entire optical depth of the hot plasma is the result 
of pair production. We have also included photoelectric 
absorption in the disc, to demonstrate the formation of the 
'reflection hump (Guilbert & Rees 1988; Lightman & White 
1988). (In our previous calculations neglecting photoabsorp- 
tion, the reflected spectra have a shape close to the shape of 
the incident spectrum at energies f <0.1. ) To calculate the 


photoabsorption cross-section we used the code of Wilms 
(1994, private communication), based on data from Verner 
et al. (1993), assuming cosmic abundances (Grevesse & 
Anders 1989; Shull 1993 for Fe and Cl) and zero ionization 
stage. The resulting average optical depth along the radius of 
the cylinder is 0.30, with the maximum gradient of the 
plasma density d r/dr ranging from 0.21 in the region closest 
to the disc to 0.37 in the inner upper part of the cylinder (Fig. 
7a, cell 2). The temperature ranges from 0.32 in the inner 
lower part (Fig. 7a, cell 1 ) to 0.44 near the upper surface of 
the cylinder. The resulting spectra for two angular ranges are 
shown in Fig. 12. The dominating slope is a = 0.81. The 
feature near e = 0.01 is a combination of the fluorescent iron 
K« line at 6.4 keV and the iron K edge, which are un- 
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Figure 1 1 . Escaping photon spectra for the case of cylinder + disc, showing etiects oi changing the distance ot the cylinder above the disc. 
Solid curve - h = 0.5; same model as in Fig. 1 0. Dotted curve - same parameters, except that h = 1.5. 



Log(e r /m e c 2 ) 

Figure 12. Escaping photon spectra for the case of cylinder + disc 
and pure thermal pair production with detailed simulation of 
photon interactions in the disc. The geometry (// = 0.5) and total 
power (/, = 50) are the same as in Fig. 10. Solid curve - total escap- 
ing spectrum in angular range cos 0>().5. Dashed curve - 
cos 0 < 0.5. Dotted curve - only the reflected component for 
cos 0X1.5. 


resolved. The deviation of the Tace-on' spectrum (solid line) 
from a power law in the soft X-ray range is probably associ- 
ated with a lower contribution due to single upscattering 
when the viewing angle coincides with the main direction of 
the UV flux. The reflection features as well as the UV bump 
are stronger in the face-on spectrum because the intensities 
of both are proportional to cos 0. 


4 CONCLUSIONS 

The comparison of our results with those of the kinetic 
equation approach demonstrates an agreement good enough 
to conclude that it provides a successful test of the LP code. 
The agreement is better than one might expect, considering 


that the different treatments of radiative transport lead to 
significant differences in escape probabilities. When we 
introduce a more realistic treatment of radiation transport 
into the LP code through the multishell representation, the 
essential differences appear mainly in the hard (f y > 1) part 
of the escaping luminosity spectrum. 

We therefore conclude that, in the case of a simple electro- 
magnetic cascade with a modest degree of Comptonization 
and moderate compactness, the escaping spectrum in the 
range e } , < 1 is not very sensitive to details of radiation trans- 
port. The kinetic equation approach as well as the PSD 
Monte Carlo scheme can be successfully applied in this case, 
especially if one chooses the correct escape probability that 
follows from our results or can be estimated from simplified 
linear Monte Carlo simulations of photon transport. 

The advantages of the LP method over other approaches 
increase rapidly with the geometric complexity of the system 
being modelled. For problems such as the accretion disc plus 
corona considered in Section 3.4.5, emission from relativistic 
jets, etc., the correct 2D or 3D treatment of the radiative 
transport becomes crucial. We believe that such problems 
can be solved with the PSD representation, but only at the 
expense of larger computer resources and with considerably 
more effort put into tailoring the code to the specific 
model. The calculations presented in this paper have all been 
made on Sun-4 and PS-486 computers, with a CPU time for 
each run of at most a few hours. In view of the statistical 
precision obtained in these calculations (as estimated from 
the fluctuations of the spectral curves shown in the figures), 
we feel that the LP method holds considerable promise for 
simulating non-linear radiating systems in astrophysics. 
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APPENDIX A: SPECIAL MODIFICATIONS OF 
THE LP METHOD 

A 1 Accumulation 

In many problems, some interactions are very soft (i.e. the 
energy transfer Ae is small) and occur too frequently to per- 
mit one to follow the evolution of an LP by taking into 
account each interaction. To treat such interactions we use 
an accumulation technique, in which a large number (F) of 
interactions is ‘accumulated’ into a single interaction. The 
energy transferred in a cumulative interaction is FAe and the 
cross-section of this interaction is <j acc = o/F. The factor F 
can be chosen so that the fractional energy change, /, of the 


incident particle is small but not negligible, say/WAf/f ~ 
0 . 1 . 

As an example, consider Compton scattering of soft 
photons on ultrarelativistic electrons. The change in electron 
energy (in m e c 2 units) is so that F~f/{e y c c ) 

and a aCL . - o T e y e c /f, \ where o r is the Thomson cross-section. 
Note that a acc is proportional to c )( , making the interaction 
rate proportional to the photon energy density. 

A 2 Time-averaged target 

There exist situations in which LP statistics are not sufficient 
to simulate a process accurately. This happens when the 
cross-section of a process is very large and the density of 
target particles is small. Examples of such processes are 
induced Compton scattering and synchrotron self-absorp- 
tion. In this case, the LP technique can be modified by incor- 
porating some features of the PSD representation. 

Let us illustrate the method using the example of synchro- 
tron self-absorption. When simulating this process, hard 
electron LPs (y c ~ 10-10 2 ) interact with soft photons. The 
number of hard electron LPs at any moment is too small to 
provide satisfactory statistics, because a hard electron loses 
energy very rapidly. This feature can, however, be turned to 
our advantage, if we store the intermediate steps in the evolu- 
tion of a hard electron LP and use the results to construct a 
time-averaged energy distribution for the electron LPs. Then 
the soft photon LP can be simulated as interacting with the 
time-averaged distribution. 

This method has the same disadvantages as the PSD 
representation - for an inhomogeneous, anisotropic system, 
a multidimensional phase-space density array averaged over 
time is required. 

A 3 Electron thermal pool 

The detailed simulation of non-relativistic quasithermal 
electrons is very time-consuming, because their energy 
changes resemble a random walk in energy space at a very 
high speed. Fortunately, we know that the result of this ran- 
dom walk will be very close to a Maxwellian distribution, 
since the Coulomb coupling among low-energy electrons is 
generally rapid enough for thermalization. This allows us to 
use a simplified method to treat the non-relativistic electrons 
as a ‘thermal pool; 

The following is a summary of how we treat this compon- 
ent. 

(1) We introduce an electron pool temperature, 7 e , and 
an upper cut-off energy, E c > kT c , At each time-step, the 
energies of all electrons in the pool are randomly sampled 
from a Maxwellian distribution with a cut-off energy E c . 

(2) Each electron crossing E c by losing energy enters the 
pool, while an electron crossing E c by gaining energy leaves 
the pool. All individual interactions of LPs in the pool are 
turned off, with the exception of pair annihilation. 

(3) The change in the pool energy A E p is calculated at 
each time-step. This includes the energy transferred in 
Compton scattering and Coulomb interactions with particles 
outside the pool and the kinetic energy of electrons leaving 
and entering the pool. The temperature of the pool is 
adjusted accordingly. 
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There is a problem associated with measuring the temper- 
ature of the pool. If we defined the temperature directly using 
energy balance and total statistical weight of pool electrons, 
we would get statistical fluctuations of the temperature (even 
including negative values) that were too large. The reason is 
that the energy flux through the thermal pool (heated by hard 


photons and electrons through Compton and Coulomb scat- 
tering, respectively, and cooled by Compton scattering on 
soft photons) is large and the total kinetic energy in the pool 
is comparatively small. To avoid large fluctuations we intro- 
duce an artificial relaxation time which forces the tempera- 
ture to adjust gradually. 




